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ABSTRACT 


The discreteness of travel-time data when using 
the Herglotz-Wiechert inversion technique results in a 
set of velocity depth profiles which will satisfy the 
data, rather than the unique model obtained when it is 
assumed that the data is continuous. The envelope of 
these models is analytically determined in this thesis 
for different types of data sets and computations of 
the velocity-depth profiles are then made for a number 
of cases. 

The effect on the envelope width of changing the 
sampling interval is determined, as is the effect of 
knowing the velocity as a function of depth in certain 
regions (i.e. in the crust or down to the lower mantle). 

In order to compare this method with other inver- 
Sion methods, it is necessary to define an averaged 
radius which is then maximised and minimised so that 
a series of extremal curves is obtained, for a number 
of averaging lengths. This type of result (i.e. the 
variation in averaged velocity at any radius versus the 
averaging length) is directly comparable with, for 
example, the resolving power versus ::adius plots obtained 
by Gilbert and Backus type inversion. In previous work 
with Herglotz-Wiechert inversion, it has not been rea- 


lised that to make any such comparison, it is necessary 
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to first average radius as a function of velocity, then 
maximise and minimise it and not vice versa. In this 
thesis these maximum and minimum averaged profiles are 
analytically determined and computations made to show 
the effect on an envelope of decreasing the averaging 
length. 

Regions of high and low velocity are discussed 
and an ‘earth-stripping' process is applied to some 
outer core data for which extremal velocity-depth pro- 


files may then be obtained. 
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Path of ray, with ray parameter p, 
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Sketch to show the effect of the ce 


low velocity channel on turning point 
of rays immediately above and below 


channel. 


A typical travel time (T-A) curve 
showing the three types of discontin- 
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BOtChR type 1° anc ii~discontinurcies in 
the p-A plane. The distance between 
the reflected and direct ray, 6éA(p) 
(as shown) is the ‘effective' range 


OL tie ~Carehn=stripping process. 


The travel-time curve obtained by the 
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Claes ce Lime 61 cere. 


BOthes this and fig. 2.60°are travel time 
data, which includes low velocity 
channels, transposed to the p-r (i.e. 
ray parameter versus turning radius) 
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top and bottom radii of the ce low 


velocity zone, respectively and the 
cross hatched areas represent the 
first of the double integrals in 
equations 2a >)" annoy (2.609% &-Figure 
2.5 shows the order of integration 


in-equation (2. 5)% 


This diagram shows the order of ail 


integration an’ equation (2.6). 


A velocity radius profile for perfect, 24 
continuous data, including a low velo- 
city zone. As can be seen, even 

under these circumstances velocity 

cannot be determined exactly as a 

function of radius below r(p'))) 
(ive the’ radius ‘at’ the top’ of the 


low velocity channel). 


The function o(p) defined in the 26 


equation for an averaged radius. 


The bounds on any model for (p,T,A) 3h 
data in the p-A plane (i.e. ray 

parameter in sec/deg versus range 

in dégs)," An. arbitrary pia)? curve 

is also shown, where, to satisfy 

the data, the curve must always be 


with the bounds. 
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Figure 


The p(A) curve, in the p-A plane, 
which will maximise ates where 
the inset refers to the change in 
endspoint 1£ q<m.. The uncons-— 
trained end point would be similar 
to that shown in the main diagram 


but with p = m. rather than eae 


sketch ofgaetamily, of gh ep CUuSVeS, 
where Hoh (oy) decreases as a func- 
tion of increasing q. Hence the 
value of m. will depend on q and 


may be constrained to u. or m.. 
J J 


A plotvinithe p-A plane of the band 
of possible p(A) curves where the 
solid lines through the boxes is the 
p(4). curve which minimises Temes 
The inset refers to the alternative 


end point value if a ee 


Radiuseveloca ty! plotlto_showethe 
effect of decreasing the sampling 
wWLothetrompebaoh tons °mandsel*®. The 
range of q values used is also 


indicated. 


Rada us=ver cea ty! plotsftory (p,T,A) data 
WLEl the, eLLect Of fixing part OF 


the model indicated by broken lines. 


A fadius=velocity, plot ywaseh] thesgen- 
velope of profiles represented by 
the solid lines and three individual 
PLOLUles, Witch are also they extre-— 


mal profiles for q=q represented 
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by broken lines. The lines across 
the envelope are the values of q, 


as shown. 


A graph of ¢(p), the averaging 
function defined by equation (2.9). 
For b<p<no, a perturbation ép 
from the mean, p = ms, will result 
in a decrease in the area under 
thee OD) Curve. UhOT ea < D< by a 
perturbation 6p from p = Me will 
result in an increase in the area 


under the $¢(p) curve. 


The values of p(A) about the point 
p = b which will, in general, 
maximise the integral, are repre- 
sented by the solid lines through 


the boxes. 


The values of p(A) about the point 
p = b which will minimise the in- 
tegral, as represented by solid 


lines through the boxes. 


The two possible ways of minimising 
the integral about the point p = b 
aSesSiOWN DY \i)jand). (2). 


A plot of extremal radius as a func- 
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TOr f= O.nSewicth Gue mode. fixed .to 
p = 7.60 sec/deg. 
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A plot of extremal radius as a 
function of the averaging window, 
(Deway re roe mW *6so°wrth the model 
fixed to p = 7.60 sec/deg. 


An alternative way of plotting 
fig. 3.14, that is by represen- 
ting the variation in extremal 
radii as a function of averaging 
length by a series of error cir- 


cles in a radius versus n plot. 


A radius-velocity plot for the 
same region as in fig. 3.8 where 
the solid line is for a zero 
averaging length. All sets of 
broken curves are for various 
averaging lengths, as shown. The 
Yelauronsco me (1.6. Ceinrtaig. 93.9) 
is shown by the solid lines across 


the curves and the values of n are 


as shown. 


These box structures in the p-A 
plane represent the bounds on 
permissible p(A) curves for in- 
version when using (p,A) data. 
The extremal p(A) curves are as 


shown peas as an arbrtrary’ curve. 


The p(A) curves which give the 
extremal values of the averaging 
integral and hence the averaged 


rads yeome(p pldsedata. 
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Figure 
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The band of permissible models, in 
the p-A plane, for (T,A) data. The 


maximising p(A) curve is as shown. 


The Minimising p(s) curve for” (T, A) 
data is represented by the solid 
line through the band structure. 
Themengd DOlnt nid), els that Lor 
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The p(A) curve which will maximise 
the integral is shown as a solid 
line through the band structure, 
with b- a the interval over which 


an average radius is taken. 


Liesminimesing: p(/\) ‘curve, for (TA) 
data for an averaging interval b- a 
is represented by the solid line 


through the band structure. 


Radius-velocity plot for “the outer 
core region, using ‘earth stripped' 
PRE alu oho Gata. 


The upper and lower limits on 1T 
for iT, p) ‘data, jshowniein sche t—p 


plane. 


Ray parameter, ,.p,,pin, sec/deg, .as a 
function, of range ,~uii, on deg... sLhis 
p(A) curve shows a triplication 
with an example of possible data, 
for unobserved triplications at 


data points, superimposed. 
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The band structure, in the p-A 
plane, for an unobserved tripli- 
Cat mon at A; - ihe cross hatched 
area shows the usual band width 
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the shaded region shows the 
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time curve, as shown in fig. 5.3 
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CHAPTER i 


INTRODUCTION 


One of *the™basie sttidies of ‘geophysics is, 
given surface observations, to determine internal 
structure. In seismology, given a set of arrival 
times of seismic waves at known distances from the 
source (i.e. a travel-time curve) is it possible to 
determine the corresponding velocity-depth profile 
in the earth? This is the problem of inverting body 
wave data, and is the subject of this thesis. 

inere™are “a °trrumber of "direct and indirect 
methods of inverting travel time data. These include: 
(1) NecglLotgs Wlecierc = ivelslOneeLreculVely inLe— 
grates a function of a travel-time (T-A) curve to any 
range and determines the velocity at corresponding 
depth (this will be discussed later in more detail). 
(2) The Monte Carlo method has been applied to T-A 
data by Press (1968, 1970 a,b), Wiggins (1968) and 
others. It involves trail and error on a large scale 
and, as such, is very expensive in terms of computing. 
The method has the further disadvantage of providing 
solutions with little real insight into the physics 
involved. 

(3) Gilbert-Backus inversion involves linear per- 
turbations of a starting velocity model until in 


agreement with known travel-time data. It has been 
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discussed by Johnson and Gilbert (1972) and others. 

Bothy (2Z)4and (3) are indiréctsmethods. That is, 
given a velocity model, time (T) and surface path 
length (A) must be computed and compared with observed 
values, to test, the model's validity. 

This thesis will be concerned with the first 
method, which is a direct inversion technique. The 
same notation as in Bullen (1963) will be used, as this 
Notation rather than\ythat of the earth flattened system 
(Gerver and Markuskevich, 1966, 1967, McMechan and 
Wiggins yell / le Wiggins, \MeMechan and ToksGz, 1972) is 
more frequently encountered in seismology (though the 
latter does have slight algebraic advantages). 

From the properties of geometric optics (i.e. for 


a ray obeying Snell's law): 
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Bro. iets Path Of Way, With bay parameter py, 
indicating various ray characteristics 
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ca is the source and receiver radii (assumed equal). 

p the ray parameter, is the distance derivative of 
the"travel=-time *curve.” It/is also thé’ reciprocal 
of the apparent velocity of a body wave across 
the surface of the earth (i.e. angular slowness). 

s(p) is the radius at which the ray bottoms. 

nN istthevratiotr/v (xy 

T(p) is the travel time for the ray with ray parameter 
Pp. 

The assumptions necessary to make equations (1.1) and 

(1.2) valid will be discussed later. 

Equation (1.1), which is*atform of Abel's inte- 
gral equation, was first investigated in relation to 
seismology independently by Herglotz (1907) and Bateman 
fL9 LO) wands lavesestmplitied by Wiechert (1910). Their 
solution, the classical inversion equation, is 


A) 
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: cosh ~ (E) aA ; (res) 


O 
1 ln — 
ase) 
This equation says that the radius s(q) at which a chosen 
ray, with ray parameter gq, turns may be determined be- 
cause the structure above s(q) has already been calcula- 
ted using the same equation. Here, A(q) is the range of 


this chosen ray. 
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This method of inversion has only recently 
been extended to include the effect of wave guides 
(Gerver and Markuskevich, 1966, 1967). Equation (1.3) 
assumesethat ta small perturbation, 6p, in the ray para- 
meter will produce correspondingly small changes, 6r, 
in the radius at which the ray bottoms and §v in the 
turning ray velocity. Obviously this will not be true 
if low velocity channels are included and the more 


general inversion equation is : 





A (q) k* x 
Xo Pie gi2 
1 ln ee | cosh (G) dA Po) | = pee 
0 k=1 a 
= atl 
tani = Ger eae (1.4) 
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(this equation will be discussed in Chapter 2) 
where each low velocity channel is assumed independent 
of the others and the summation represents the effect 
GE vChaniicl sale torougiek*. JAll other isymbols jare, as 
previously defined, where the subscript k refers to the 
ka low velocity channel, as illustrated in fig. 1.2. 

| A number of assumptions must be made in this 
inversion; 


(1) Both source and receiver radii are equal to ro: 
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Fig. 1.2. Sketch to show the effect of the ane low 


velocity channel on turning point of rays 


immediately above and below channel. 
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(e2) Only one ray type and velocity is being investi- 
gated. 
(3) The velocity v(r) is non-zero and may be discon- 


Pinole -aseCahethe sti rst and second derivatives ,.but 
only at a finite number of points. Elsewhere the 

second derivative must be continuous. 

(4) It is assumed the complete T-A curve is known 
exact lyveui.e. gall values ot sp girom gq, £o ri/Vo Fay) 

are known). 

(5) Pneresi Ss novambiguitypofts2n7 inh ,.except ipossi- 
bly in the inner core. This is reasonable as rays which 
have path lengths greater than 27 will have much reduced 
amplitudes. 

If all the above assumptions are true, then 
knowing p = Wotatei ws Op iteisepossible  toecalculate 
velocities near the surface radius ro: As p decreases, 
that is as deeper rays are sampled, corresponding velo- 
cities may be determined because the velocity structure 
above tle ray being sampled has already been found. 

Practically, data will be discrete rather than 
continuous, so that the T-A curve will be represented 
bysavseG of points rather than a continuous dine. “That 
Lisi es moving from one range to the next there is an 
ambiguity in the velocity depth structure "seen" by 


rays at these ranges. For the Herglotz-Wiechert 
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inversion method, this discreteness of the data means 
that for a chosen ray (with ray parameter q), there 
will be an indefinite number of velocity models above 
the radius (S(q)) where the ray turns, which fit the 
data (because of the above mentioned ambiguity). 
Hopet utly“toreach' chosen ray, St’ will be’ possibile’ to 
select the velocity models which will maximise and 
minimise s(q) (and hence maximise and minimise the ray's 
Vevrocyty sat tite tuning port * tow ‘ve Foeci#ty-"ones are 
temporarily ignored). This thesis is concerned with this 
PLovlen. 

As will be shown in Chapter 3, discrete data in 
a I=Ccurve’ Corresponds to-anvenvelope am the™p-A“plane 
(o =—"d't/ dn) and it is’ possibile’ to“analytically *determine 
Pu curves 1 tliils plane witci: 
(1) give extremal values of the integral of equation 


(1.3) and hence extremal velocity-depth profiles, 


(2) lie completely within the p-A envelope, 
(3) satisfy the travel-time constraint 
A(p) 


The presence of low velocity channels presents 
GitfireuLties and though at is still possible to, Lind “the 


extremal paths of integration, the complete analytic 
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solution has not yet been obtained. 

Regions of high velocity gradient have been 
treated by McMechan (1971) and will be discussed in 
Chapter 4. 

Arbitrary p(A) curves in the p-A plane corres- 
pond to non-extremal velocity-depth profiles which 
will certainly touch the extremal curves but will 
drop away again (this will be shown in Chapter 3). 

For example, at a particular radius, an arbitrary 

model will actually define the maximum velocity. 
Bessonova-et al (1970) and Keilis-Borok (1970) define 

an average radius which is effectively the turning 
PactUcseeLolrciimavelaged. VelOClucy proLile. (ne COLLes= 
ponding average velocity-depth curve may also be maximised 
and minimised, as will be shown in Chapter 3. 

A trade off curve may then be computed for the 
width of the averaged velocity-depth envelope versus 
the range of p over which the averaging is taken. This 
will provide a direct comparison with the Gilbert-Backus 
results (which is intimately concerned with the trade- 
off of resolution versus accuracy) and the Monte Carlo 
method (which will give a low "density" of models through 


extremal points). 
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CHAPTER 2 


BASIC METHOD OF INVERSION 


Initially perfect, continuous travel-time data 
will be considered, so that certain features of the 
travel-time curve may be determined. Regions of high 
or low velocity gradient introduce various anomalies 
into travel-times, as illustrated in fig. 2.1. These 
anomalies are: 

(EL) Am@cCoLbiicacton, whici=iaea region of high velo- 
eCity graadient. “For continwioustdatay, «triplications cause 
no problem in inversion and have been investigated by 


McMechan (1971) for discrete data. 


(2) eae yper le discontinuity, which is due to a dis- 


continuous decrease in velocity gradient at radius y (K) 
such that 
nie +O) < Hee = .0) 
(k) ans) a. 0) 


where (4r 4 0) is the radius just above and (r 


ele An example of this structure is 


that just below r 
the core. 


(3:) Ajtype Il discontinuity, as with the low velo- 


city zone, is caused by a continuous decrease in velocity 


gradient at Che such that 
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In practice (2) may not be a clear discontinuity 
Tomiw because tof siicaring by diffraction. A third type 
STeOVecOntinuity, sl. ,emist also be considered. This 


iseaycombanation of types 1 and IT, such that 
ie()Pes) ree) (2-82 ) 


PUG ee eu eevee ue) 
This type of discontinuity may be a result of uncertain- 
ties in types I and II and it has been suggested (Bolt, 
1972) for the D" branch, at the core-mantle boundary. 

oy Demi mOG el OLrocontinuity may. be inverted 
MOLCmeaS Ly etieanra Voce sll discontinuity and it is 
therefore important to determine their differentiating 
characteristics. Before proceeding further, it should 
be noted that the parameter t, defined in fig. 2.1, may 
also be written in integral form 
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Thus, it is very difficult to determine the type of 
discontinuity directly from the travel-time curve and 
amplitude properties must be inspected instead. 

Using wave theory it is possible to show that 
for type I, the amplitude of a diffracted signal decays 


exponentially for’ A > Noa 


with a cube root dependence on 
frequency (equation (A1.3)). Similarly, for type II, 
the diffracted signal amplitude decays exponentially 
for a> al) but is independent of frequency (equation 
(A1l.5), j=0). This result agrees with that obtained by 
geometric considerations (equation (A2.3)). Fora 
type I' discontinuity we find that the amplitude beha- 
viour is the same as for type II but with a decay rate 
which depends on the boundary conditions (equation 
0) J 

It is assumed that the velocity model to just 
above the low velocity zone has been inverted and is 
known. Ideally, it is possible to determine the type 
of discontinuity through amplitude inspection, though 
in practice amplitude data is difficult to analyse. 


(k) (k) 


For both type .' and -diydiscontimuities. 4 and T 


are infinite and obviously the whole curve A(p), as 


oe cannot be observed. However, it is 


(k) 


p tends to p 


possible to extrapolate the model to r , uSing equa- 


(k) 


tion (A2.1) and hence calculate T (equation (2.4)). 
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Similarly, using equation (A2.2), A(p) can be extrapo- 
ated, Lor p tending to alto from above and below and 
hence an estimate of aou obtained. This then provides 
a value of sae) where 

tk) 

2 3% 
a Ee ia (nF oN ar 

a) 

A second constraint, 5A Sle) , which is illustrated in 


fig. 2.2, 15 indeterminate for a typefll discontmuity 


(k) 


but remains finite as p tends to p for both types 1 


enael's or the Wattem, the wanges themselves become 
intinvte about the discontinuity though their difference, 
vie 


remains finite. This parameter is not directly 


observed but may be calculated as 
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ae one Ds 


o je) teks 
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These measurements 6T and 6A , whether directly 
observable (type I) or estimated by extrapolating of 
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pe gl lustratedin fig. 2.4, Which shows the "travel 


al 


-oqge13x9 ed A695 (qg)A .(S.SA) mottsupe paias yirelimia 
bus woled bas svode moxt re ot pnaibases g toz ,betsi 


asbivorg noid erit benistde [5 30 sismises us soned 


‘exedw (0) v4 Yo ovisv 5 
, OD, 
46 € Sttg = *nkS= se ag 
(#5 
oi betdsageulli ai n>idw tay a8 ,tnisitenods Draovsa A 


ysiunisnessiS IT sqyji 5s 102 sasfhimxstebni ei, S/S .ees 
I esqys dAfod tot A) 6, o3 ebnet q es etinet enismet Jad 


emoosd asvisamedt eopnsi edt ,1tettel sdt 107 .'T Bas 
} = 


,S0nesstILb rziocds Apu yiiunisnoo2zib sit tuods stialtat’ 


yisoeirb toi ei zetemexsq etdT .stini? anremst ins 


26 betsivolso ed ysm stud BeyazsRdo 


ve 
: 0 


= 
7 ‘oa ays (A) S 


bits 
(gy nae “ne 


Yisos1ib rsitedw ne bas ‘“) 5, exnsnsiuesom oeedtt 


20 








Pigs 


ere er 


15 





Pitustracion OL the -"eartii-stripping 
CONCeEDL TOr a ray eltering a type If 
row’ Vvelocrty regione in Chis’ Case “tie 
rays used are PcP and PKP (Macelwane 
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Eimescurve for such a “stripped” earth. An “earth-— 
stripping" process (Macelwane, 1951) removes the 
reflection time, which effectively removes the layers 
abovemtne discontinuity, aS shown in fig. 2-3 in the 
D=Agplane. § Thus, a Te ee ane? (p) curve is obtained 


which will have no observed values below REALE 


Pee Perigo Satie cres iunay be possible to 
fill the lower part of this curve with other data, 

for example SKS for the outer core. If this supple- 
mentary data does not exist, the situation becomes 

a special type of discrete data problem, which will be 
discussed in Chapter 3. Witthethais—technigque, it is 

most convenient to assume that there are no triplications 
om low velocity zones” within thesregron=of=the tow velo- 
city zone (i.e. velocity increases smoothly with depth 


(k) ne 


between r and r Multiple zones could be inclu- 
ded but will not be studied here. 

The above process is not possible for a type II 
discontinuity because Ny Tey is undefined. Instead 
the results ob Gerver and Markuskevich (1966, 1967)7 
which generalise the classical Herglotz-Wiechert inte- 
gral to include low-velocity channels (equation (1.4)), 
must be applied. These results will now, therefore, be 


derived. It should be noted that to compare the follow- 


ing with the results of Gerver and Markuskevich, the 
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Pp reflection from 
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reflection from 
r (2) 
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Pidwacese  BbOth type 1 sand Lisdiscontimuruies: Tene 
p-A plane. The distance between the reflected 
andsdirect. rays, sA(p) (as shown) as the 
‘effective’ range of the “earth—-stripping~™ 


process. 
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Fig. 2.4. The travel-time curve obtained by the "“earth- 


stripping" process, as illustrated in fig. 2.3. 
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transformations from earth flattened to spherical co- 


ordination must be used. The operation 


Oo 

ce Ree ee 
| seer 
g Ped 7) 


1s applied to A(p) (equation Al.3d)) := 
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The right hand side of equation (2.5) is equivalent to 

an integral over the area p £9) (aq) ingthe, p-teaplane, 

as shown in fig. 2.5, where the total area between the 

DOLNCS = (Sia), 1 (e)) and (ror q) is defined to be D(q). 
k 


k=1 
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Both this and fig. 2.6 are travel time data, 
which includes low velocity channels, trans- 
posed to the p-r (i.e. ray parameter versus 
turning radius) plane. In both figures 

Tye ae are the top and bottom radii of the «eh 
low velocity zone, respectively and the cross 
hatched areas represent the first of the double 
inteqralsisin equations. (2.5) cand. (24.6)... (Frag. 


2.5 shows the order of integration in equation 
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Fig. 2.6. This diagram shows the order of integration 
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In general, 


oS ote UNS eat 5) Lena # 
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and hence the first term on the right in (2.6) reduces 
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integrating the left hand side by parts gives 
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and so the equation becomes 
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Obviously, for pit? saq Soe “tithe” classical mesuitras 
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obtained, where the summation term is ignored, and 
for "perfect" data this will produce a unique velocity- 
depth profile. As explained above, for q tending to 


pit), the upper depth, (x {1)) at the low velocity channel 
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is*known, but ts “extent is not; asi it depends on how 
the velocity function varies within the channel. 
thusyeeven for "perfect! data, it is only possible to 


impose limits on the velocity model for gq < pane 


(1) 


using the constraint 6T Cog.) 82iet7-) eel SewL.Lbe 
discussed in more detail in Chapter 3. 

With discrete data, p(A) is not known exactly 
and hence Tet) cannot be evaluated explicitly. 
However, if maximum and minimum values of To) (aq) 
can be found, while still satisfying the data, this 
will give minimum and maximum velocity and radius 
values and define the envelope mentioned previously. 
The problem of inverting low velocity channels of type 
II can similarly be specified as finding the extremal 
values of Th (eq) which fit gehe data Kcontinuous or 
discrete). 

When considering low velocity models, numerical 
instabilities are encountered, in that thin, very deep 
low velocity zones cause only small changes in T(A) 
(Kewlis — Borok, 19/0, Bessonova ct al, 19/0) Hence, 
rather, than .calculatang- the struc. proltiie “allwaverage 
structure may be used, which will be more stable. In 


radial variables, the definition equivalent to Bessonova 


et al's (1970) is 
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r(p*”) r 


AO VNGTOC Eyer usiprorile £or perrect, con 
Siuuous date, ancluding a iow velocity zome. 
As can be seen, even under these circumstances 
velocity cannot be,determined exactly as a 


(1) 


function of radius below r(p ) (1 se. the 


radius at the top of the low velocity channel). 
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ilial R(a,b) = pen | in s(n) dn 
a 

where a < p s b is the averaging intervaland, as defined 


earlier, s(n) is the radius at which a ray, with ray 


parameter q, turns. We recall that (equation (2.7)), 


Ne 
in sq) = ae) - ya ieee wise ; 
k=1 


ce) 
and thus 
: k 
Dale pea | (Po @ity 2 1° tated 
O s k=1 


The first term on the right hand side becomes, on 
substitution for 1 ') (g) and reversing the order of 


integration == 


n b 
oO 
| A(p) {cos ? = - cos? >} 4p = | hp) cos a dp 
b . a 
he 
= | A(p)o¢(p,a,b) dp 
a 
where 
Sloe “lb 
— = COs = bus < 
cos = c = p ne 
¢ (p,a,b) = 
cos’? 5 Bh is) bs) 


(as shown in ifig. 2.8). 
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The function $¢(p) defined in the equation 


for an averaged radius. 
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Hence 


egies 
A(p)o(p,a,b)dp + | eta (q) dq]. 
kaw 


Bessonova et al (1970) introduced this average to avoid 
the numerical instabilities of solutions including low 
velocity channels. However, it will be used here in 

the interpretation of discrete data. This problem will 
be looked at in the next chapter. The envelope of all 
models defined by isco) and me (ay will be discussed, 
as will the maximum and minimum average models defined 


by R 


max (2rP) and Rain ‘@rP) - The former is actually a 


special case of the latter with a = b, where (a,b) may 


be thought of as a resolving window. 
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CHAPTER 3 


DISCRETE DATA 


inereality; rathers than an infinitesamount of 
data, there will only be a finite number of data points, 
assumed #60 be, part fof a “continuous function. (In this 
Chapter, it will be shown that even though the condition 
of knowing a complete T-A curve is not met, it is still 
possible to use Herglotz-Wiechert inversion techniques 
and obtain a set of velocity-depth profiles which satisfy 
the data. 

First, extremal velocity-depth profiles will be 
calculated for the forward branch of T-A curves only, 
where dp/dA < 0. Regions of low velocity gradient will 
be discussed in the next chapter and reversal regions, 
antler have been considered by McMechan and Wiggins (1972), 
will be reviewed in Chapter 5. 

The most highly constrained case to consider is 
when phase velocity (p) as well as travel time data is 
known. That is ee T(A;), ee p(d;) are exactly known 
for certain stations Ay i=l,..., where it is assumed 
that Ty = 0, AD ea 0 and data are arranged in increasing 
Order: 


Though information is only obtained for a finite 


number of points, the condition that the gradient of 
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T = T(A) must be a monotonically decreasing function 
of range provides a further restriction. That is, for 


the linear interpolation gradient m, m <oeNlea 


itl Lf ia Y 
where 


For (p,T,A) data, the upper and lower limits on p may 
be further reduced, as ms represents an averaged value 
Of ray Parameter p between (T;,A;) and (Tey Agar? 


while u; and U;,, are the p valuesat iy. and 


(Tor Aayy) respectively and therefore, as dp/dA < 0, 


oe eee Se 5 
a0 1 it+l 


The constraints on the data may be summarised 


as 
(3) a ae fOr A; a 
(2) p= u; 1m ele A = A; 
dp 
(3) ai 0 
A 
(4) | p dAv= 1. 
0 


Equation (4) may then be modified, to 
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Constraints, (1) and (2) state that in the p-A plane 
any valid p model must lie within independent box 
structures as ilitetrated in fig. 3.1 which indicates 
an arbitrary p curve and™the bounds on any p model in 
thnegp-4 plane. 

The problem then becomes that of determining the 


(Oo) 


maximum and minimum I q) (as defined in equation Zo) 
subject to the above constraints. Obviously, the 
simplest way to maximise the integral is with p = U;- 
However condition (4) would not then be valid and in 
fact will only be met if p is evenly distributed about 


the mean, mM, 


A 
GO ey | ae e an 
0 


and 


Thus any small perturbation of p, ép, from p = m, will 
result in a decrease in the total area of the interval 


aS indicated by £19. 3.2-. 00 £Orep = Mm: the integral 
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the bounds -Oneany models Lor (p,i,A) data 
Diecierp— eo lancuil.c. ba yvepalanle Comal 
sec/deg versus range in degs). An arbitrary 
p(A) curve is also shown, where, to satisfy 
the data, the curve must always be with the 


bounds. 
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larger of the two increments 


vA 





Bigs 3.2. Amsketch of cosh + py qvas a function vot yp) 
showing that any small perturbation, op, 
away from the mean Ms» will result ina 


decrease in the value of i?) (aie 
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will be the global maximum as shown in fig. 3.3. How- 
ever, the end points must be considered separately. 
The end point A(0) is quite straightforward as p = my 
will produce the maximum, as in the general interval. 

At A(q), the integral will be maximised by taking 
p= m in the interval containing p = q, and, beyond 
Aig), by taking=psequal~to-its.minimising value. The 


relationship between A(q) and m. ZS E£ound by the tollow= 


ing equation, which satisfies condition (4), 
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Hemecany be ceen fbom tag. 3.4, the value of m. will be 


dependent on the value of q, and as A. SIN ep) SN 


ms is constrained above and below to a, and Mh respec-— 
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or q 


The p(A) curve,in the p-A plane, which will 
maximise mo) 7G) where the inset refers to 
thes change simp ends pOlncetesds— ae The 
unconstrained end point would be similar 

to that shown in the main diagram but with 


p=m, rather than tee 
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or q 


The p(A) curve,in the p-A plane, which will 
maximise Ma) where the inset refers to 
the change in end point if q < Hs The 
unconstrained end point would be similar 

to that shown in the main diagram but with 


p=m, rather than eer 
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Pig. 224. Sketch of a famuly.on t £9) (q) curves, where 
tO) (gq) decreases as a function of increasing 
gq. Hence the value of ms will depend on q 


and may be constrained to ue or ae 
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For een <0 , the constrained maximum is at 
cle 
Dea ei andstors (ol. /om vase > 0, the constyained 
J J m.=u. 
41° eel 
Maximum is at p = u.. 


To maximise Tae it is necessary to perturb 


p from the mean by as much as possible, say to 
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as shown in fig. 3.5, a p-A plot of the band of possible 


models where the solid p(A) curve through the bands is 
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Figee 3). “A plot =dn*the p-A plahe of the band of possible 
p(A) curves where the solid line Ehrougim wie 
boxes is the p(A) curve which minimises eo en. 
The inset refers to the alternative end point 
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that which gives the minimum integral. 
The end points must again be considered separately. 


For A(0). there. are two! cases: 


Gis) ee (i.e. no surface velocity constraints) 
then 
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ey ae -1 Y2 
I = (65 Ay) cosh q " 


(ii) uy yo 
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At A(q), the global minimum will be 
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Thus, for ue <aq< Bsc the integral - must be cons- 
trained to zerg, as *shownsinm figs 3.5. 

Using these conditions on maximisation and mini- 
musation a program was written to find extremal velocity 
dépthemodels for given (p;,A4) data. Test data (from 
the Herrin tables) was used to determine a number of 
effects. First, the sampling interval was varied between 
approximately 1, 3 and 5 degrees for this standard data 
andwextremal profiles were calculated for 9.0°= gps; 10.0 
qen0elecec/deq anternvalseor gq (this corresponds to a 
radius range of about 5600 to 5800 km). As can be seen 
from fig. 3.6, decreasing the sampling interval decreases 
the width of the velocity-radius envelope. This is as 
expected as in the limit of continuous data the velocity 
depth profile would be unique. 

Tests were also made on the effects of errors in 
T and A measurements but the results of this were in- 
conclusive and haye not therefore been presented. Any 
meaningful peewee rom this type of experiment would 
require systematic variation of both T, A and p measure- 
ments (both dependent and independent variations) over 
different interval lengths and at different points in 
Ge data Sec. 

The model was also fixed to different depths and 
the resulting profiles compared. That is mikigohen cee) 


assumptions were made about the earth's internal structure, 
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enly that 1ts surface radius was 6371 km, Mand using 

the standard data the velocity-depth envelope for the 
crust and mantle was computed, as shown in fig. 3.7. 

It was then assumed that velocity was known as a func- 
Giron Of radius jm the crust so that co a depth of about 
SO0ekm (1.6. madius at 6322 km) the velocity-depth 
profile from the inversion procedure is unique. Below 
that, a series of profiles is again possible and an 
envelope of profiles is obtained. Thirdly, the model 

is fixed to a depthaot approximately 850 kmi(a8e. it#is 
assumed that the stmucture of the upper mantilesis also 
known). The last two cases have been superimposed on 

the radius-velocity diagram of the first case, as is 
shown in fig.3.7. It should be noted from the diagram 
that the effect of fixing any part of the profile only 
reduces the width of the envelope in the immediate vicinity 
of the known part of the profile. Thus even knowing the 
structure of the upper mantle the band of possible velo- 
cities for given radius close to the bottom of the mantle 
remains unaffected. This is also the conclusion drawn 

by Wiggins et al (1973) though their method of Gnversion 
is different and uses few time constraints. However 

the data used here is meant only to illustrate the theory 
and may not, therefore, be compared too eritically wath 


the results of Wiggins et al. 
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It is also of interest to note the behaviour of 
individual profiles as they move away from the extremal 
curves. It should be remembered that a function of 
p(A) which will maximise the integral, ree) (aye for one 


value of q, say ce Wlll not do so 2or other GQ values 


x 


(i.e. for other rays). In fact the way in which the 
v(x) profile moves away from the envelope edge (a7 bili 
case the minimum curve) will depend on the position of 


gq with respect to the p parameters at the end point. 


max 


* 
F i fe Pe 3 7 
Or example, if at the end point Wu u. tO ay 


(and p = m5; for all other intervals) then for 


< = i 
Cpreenge sd. i ery these values of p will also be the 


maximum ones so that the model will be on the envelope 


LOL). SeCe s ih The system is also oscillatory as 


max 


the same conditions will be met (i.e. the given values 


Gisp will be those providing the maximum integral) LOL 


x ; "Ok : fp 
Us4, 2 48 m; foreaii@i <9 |. soimilarly,. 10 fOr. q— Gye 
m: = he then this set of p values will be the maximal 

g < < <a: gacsenr tes 
sev for eee Sea OS we and for U4, 742 Me ns 5 


as m., < ms < las then this function of p(A) will only 

give the maximum integral for u;,, ore Ges mss pitanja 
Obviously then, the individual characteristics 

of the v(r) curves will depend on the distripution of 


data with respect to the values of q taken. Computa- 


tions were made with the same sample data as used before, 





cb 


to woiveded sit” ston o3 teoxstnt to coals ef JT 
ismettxs sdt+ mort yews svom yes as esltiosq Isubivibat 
to aottonyi s 35cd¢ bexsdmemex od bivota JI. .eevigo 
eno 1t0% vbpy 80s .lsxpeotni ot seinixem Iliw dotdw (Ajg 
eoul sv p isn%o 1o> 02 ob gon Lliw xom= YS® +P io ovisv 
oif3 noinw nc ysw oft gost al. (eysz 1930 x02 .9.4) 
eif_a ni) spbs sgolevas oft moxt ysws asvom Siizouq (a)v 
to motsiecg sdi no brheqsb Iliw (svivo mUminim sae Sago 
.Jmioqg bao of3 35 2tstemsisq q sit o¢ tosqasr dtiw cau” 
xem? ~P 102 Pad = i jniog bns edt 36 21 ,slamexe FOF 
to2 AD ACRE SITES x9fito IIs rot pe =a oan 

efit sd oals iit q to asulsv sears [- ~~ zpe2 steal 
egofsvas efit no sd ‘LEBe isbom edt 3st oe 29no0 momixem 
as yiossllioao o2ls et meteye sdT “parm 2pe xem? 702 
soulsv asvip oft .9.t).t9m sd iLfw enofsibaoo smee ong 
107 (fexpstai mumixsit sit pnibivorg seeds ad ae q to 
‘xem? =p x02 Tf ,yvixslimie .f > i Its toi = 
famixem sit sd Iliw esalev q to toe. aidéoneyis ps * 
a Et eee Sa 
Yltto ILiw (A)q Xo not torus? aint nos Ere > ‘m > tt 21 


Fane ¢ sym 2p 2 ase xox iszpeait: muni era art - 





- ps oe 
a 





44 


with the profiles fixed to a depth of 1325 km. In 
Cm. 6 Only alcomall section Of the total profile 


om Plustrated (i. @.9.6..0 “<s0 6.9)8 so. Chat the 


< 
naxx. = 


detailed characteristics may be observed. These 
individual curves do show the above properties. For 


example, at 2 Bre Ont Mm, < m: < u. and the curve 


moves away from the minimum quite fast but returns at 


Genesee 6. 6 MWh Che corresponds tO (u. 


< < Gi, we 
oa giall VP de 3° i 


ak je etre Seay es 

The individual curves of fig. 3.8 are possible 
velocity models which, though they do touch the 
extremal curves at certain points, in general lie 
within the envelope. The extremal curves (which are 
not themselves possible models) represent the largest 
range of velocities at any one depth which the given 
data will allow. Hence, it is unreasonable to compare 
these extremal envelopes with results obtained by 
Gilbert and Backus inversion (Johnson and Gilbert 19 72)% 
as has been done by Wiggins et al (1973), because the 
envelope from this method of inversion will obviously 
be broader than the averaged uncertainties from other 
technigues. In order to compare uncertainties in this 
method with those of other inversion methods, it is 


necessary to use an averaged radius, as defined by 


Beesonova et al (1970). 
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This averaged radius, given by equation (2.8), 


Ls 





ln R(p,a,b) = ln s(q) dq 


b-a 


he 
@ -— 0 


where averaging is from p = b to p = a. Obviously for 
b = a, the averaged radius is R(b,b) and the extremal 
Fada. Rnax (Prb) and Rain ‘PP? are as discussed above 


(i.e. equal S max (P? and Sai 


By (Pde) For the case under 


Gonsideration, that is no low velocity zone existence, 
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ce iE ish -l b 
Coss) =" cos’ = ley < 
Pp Pp = "o 
o(p,a,b) 
Sone. a= (Dp) <7 b 
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To find the extrema of R(p,a,b) it is necessary to 
consider the extrema of $(p,a,b)A(p). 

In general,ydp/dA < 0 is required, so that any 
perturbation, 6p, away from the mean, p = Ms Wall 


decrease the total area under the A(p) curve for that 
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iceanrvals se Simitarlyyetear b <p < ihn AG /dpr= 0" and, 
aspcan be seen from fig. 3.9, a perturbation, op, 

from the mean value of p (i.e. m;) will also decrease 
the area under the ¢$(p) curve for the i interval. 
Therecore, for b p< nh, p.= m; will maximise ¢(p) A(p) 
in the eae interval. However, for a < p < b, do/dp > 0 
soithat, as’ shown in fig. 3.9, a perturbation dp from 
p= ms will ee thesarea under the ¢(p) curve for 
that interval. This confl:.cts with the effect on A(p) 
Of perturbing p. For the data iconsidered, the possible 
Varia Sron=in=et-(p)==GresS irom A. s6(®) Asy? is smaller than 
that in $¢(p), which varies quite rapidly, so that the 
maximum possible perturbation of p from p = ms which 
satisfies the travel time constraints, will maximise the 


wneeqre tor a < pi < b>. In the case of (Dp, Ly) data, 


this means that 
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larger increment 





Fig. 3.9. A graph of $(p), the averaging function defined 
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Conversely, tor b Gap < Nor the largest possible 
perturbation of ép from the mean will minimise d(p)A(p) 


So that the minimising values of p are 


(oe ih Go ke 6s. 
i 4 


harps | O 3 — Av4 
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Forya < p =< b¥ the effect on ${p) will predominate so 
that p = Me will minimise the integral in the yee 
interval. 

As before, the end points must be considered 
separately. For both cases, MOyEAsMtreated@as Gor 


(q) and s in (D> liga Usa S be < u; and a > u; 


Smax m +1 


then A(a) is also the same as the end point At qj ban 

the previous discussions. The most general maximisation 
of the integral is as shown in fig. 3.10 where A(b) 
Corresponds to A(q) for bo< pe< 1 andy (0) itor ta < pi. 
However, this contribution may be further increased by 
decreasing m; and) 1ncreasing e. until in the limit 
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* 
m, may be determined numerically by substituting 0: from 
(3.1) into equation. (3.2) and then using the condition 
* * * 
oT, /dm, = 0. Both O5 and m; are subject to further 


constraints: 
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Ci) be ms u; 
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(2) A.< 0; s Me 


The right hand equality of constraint (2) sets ms = ™,, 
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The general condition for a minimum in the inter= 
val Containing b is as shown in fig. 3.11 where Q; and 


* A 
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The: valuessor p(A) about the ypointlp— b 
which will minimise the integral, as 
represented by solid lines through the 


boxes. 
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The absolute minimum for this interval can be numerically 


determined. The two limiting cases are:- 


(1) us m* m. 0. A; 
(2) lq = Ae m+ = u; 6 = 05 , 


th interval IL 


liinie et Usays ie walbivonky atiect thert 
Mikeecero ee) Len m. = a rather than m.. 
a a i a 
Thus the function $(p,a,b)A(p) may be suitably 
maximised and minimised to obtain the extremal average 
radii. Using the same model as for the calculations 
of the maximum and minimum profiles (i.e. Sed) and 


Ss (q)) fixed to the same depth maximum and minimum 


min 
averaged radii (Riax (41) and Rin 2")? were computed 
for a number of window lengths at fixed n. The values 
of n were chosen so that the same section of the profile 
was observed as for the non-extremal profiles of fig. 
Jeee ise. Gels 75 6.5) and the window length, b — a; 
was varied from 0.1 to 0.9 sec/deg by 0.1 intervals for 
each value of n. The values of Rnax (arb) and Rain 627?) 


= sui j f athe 
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window length was not increased further as at a length 
of 1.0 sec/deg the averaging is over about 200 km and 
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large scale variations and the values obtained would 
no longer be meaningful. 

It should further be noted that the averaging 
is a logarithmic averaging over a non-linear curve so 


that the separation of R (a,b)) and Rk 
max m 


in (2 1P) widsl not 


be a linear function of b-a. For example, if ona 
very steep section of the velocity curve, then the 
average will be weighted t) one end of the averaging 
interval so that the values of the averaged radii may 
be higher than the non-averaged values. However Li tchis 
is the case both maximum and minimum would be weighted 
in the same way and both should move together. 

As with the non-extremal profiles, the position 
of n with respect to p (and hence of a and b with#respect 
to p) determines the behaviour of the averaged radii in 
moving away from Lop MY) and Sine so that 22 an 
the non-averaged case a particular profile moves slowly 
away from the extreme then it is to be expected that in 
the averaged case, the radii will move slowly away from 
their values at zero window length (i.e. jelghehe. (eee ke 
separation decreases slowly as a function of increased 
window length). This phenomenon may be observed in figs. 
3.13 and 3.14 as compared to fig. 3.8. Thus for q=6.5 
in fig. 3.8, the s(q) curve moves from the extremal 


value slowly and, as can be seen from fig. 3.13 (which 
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is a plot of averaged radius R versus window length, 
b-a for n = 6.5) the separation of the averaged radii 
also decreases slowly with increased averaging length. 
However, Rene ="6.3 in tig. 3.8, the non-extrnemal model 
moves away from the envelope quite quickly and the 
separation of the averaged radii also decreases quite 
fast with respect to averaging length; as can be seen 
from fig. 3.14 which is again a plot of averaged radius 
R versus averaging length b-a but for n = 6.3. Though 
the rate of decrease of envelope width varies for both 
n = 6.5 and n = 6.3 the width of envelope decreases by 
approximately 40%. 

Fig. 3.14 may be plotted in another way as in 
figu. 5.15 whichois-a ploG’ oF averaged radius against n. 
Thus for n = 6.3, the variation in envelope width as a 
function of averaging length may be represented by a 
series of error ellipses. Obviously there will be an 
optimum value of b- a (and hence Rmax (arb) and Rain (212) )s 
as for large values of b- a, -though~the separation 
between extremal radii is small, the area of the error 
ellipse is larger than at some intermediate values. The 
envelope for a = b has been superimposed on fig. 3.14 
and it can be seen that those cases of larger b-a 


values do not follow the trend of this envelope. 
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Yet another way of expressing this data is in 
Seroius- Velocity pDlOtpas anh fig. 3,16. Here the 


averaged velocity has been defined as 
Vilay oy =) 2Relan be) 7 wart Db) 


The region considered is the same as that of fig. 3.8 

and the outer envelope is for b-a = 0 (i.e. eee 1d) 

and Simin ‘a? > The series of curves inside this envelope 
represent various values of (b- a), as indicated on the 
diagram. As previously indicated ,@it can be seen that 
the effect of increasing the window size decreases 
especially beyond b-a = 0.4 sec/deg and that the changes 
in envelope width are not uniform with regard to different 
values of n. However,cat b-a = 0.4, the averaged 
envelope has been reduced to 60% of the initial envelope 
for all, values offm considered. 

As was discussed earlier, comparisons with other 
inversion techniques are only valid if an averaged en- 
velope is used, rather than the non-averaged extremal 
curves described at the beginning of this chapter. This 
type of averaging has not previously been done when 
using Herglotz-Wiechert inversion techniques. In par- 
ticular, the reason Wiggins et al (1973) were unable to 
obtain results consistent with those of Johnson and 


Gilbert (1972) is because they were trying to compare 
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two different types of results. Averaging extremal 
curves (as oppose to maximising and minimising averaged 
curves) will smooth out the envelope but will not reduce 
its general width. The results of Wiggins et al 
represent the most extreme velocity variations at any 
one depth, for the given data, while the resolving power 
versus radius plots of Johnson and Gilbert (1972) are 
plots of variation in averaged velocity at any one depth. 

For this reason, the type of results shown in 
fig. lopare important, )in, that they provide the only 
way of testing whether this method is as good as other 
inversion methods. However at this stage, it is impos- 
sible to do a quantitative comparison with other tech- 
niques as no real data has been used as yet in this 
inversion method. 

The method of maximising and minimising both 
velocity-depth models and averaged models discussed 
so far in this chapter may be extended to other types 
of data. As the arguments are basically the same, 
they will only be outlined in the following. 

The second type of data set to be discussed is 
(Mec tta.. = Liate Ls ipa p(A;) for a finite number 
of stations Avs ley ee AS eG tLe sud cana availa- 
ble, the constraints on any p(A) curve are reduced 


because the travel-time restriction (4) is no longer 
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applicable. Thus the p(A) curve which will maximise the 
interval may be taken along the upper edge of the box 


structure in the p-A plane (as shown in fig. 3.17). 


Theres 


Similarly the path along the lower edge of the box struc- 
ture (fig. 3.17) will provide the minimum integral where 


i a Ls Asal 


For both cases end points need not be considered separate- 
ly. 

Maximising and minimising the average radius 
R(a,b) is also considerably simplified, as with no time 
constraints it is only necessary to maximise and minimise 
A regardless of the shape of >. Hence the maximum and 


minimum averaged curves will be obtained from 
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_mMaximum curve 


_-Mminimum curve 


arbitrary curve 


BELG. aseLjeyebnese box structures Lipchcap =) Blane 
represent the bounds on permissible p(A) 
curves for inversion when using (p,A) data. 
The extremal p(A) curves are as shown , was 
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respectively, as shown in the plot of p versus A in 
figao,le. Thus this type of data would be easy to 
program though it has not been done in this thesis. 
Finally, sets of (T,A) data must be considered. 
Thag is, T; = T(A;) for a finite Avy Poly e e e re 
no p data is available and so the restrictions in the 


p-A plane become 


ai 4 Ss S Ul. : 
cell) We Does ae, for PE Se aay 
(2) ps us for A= A. 
dp 
33) ane 
Hep 
(4) | pacar. 
A, 
a) 
where u. = mM. = wW. so that rather than independent 
a i=l tha 


box structures in the p-A plane, there will be a band 
of permissible models, as shown in fig ese. EL 
addition, there may be some phase velocity measurements 
which will further restrict p in the intervals to which 
they refer, though the limit when p measurements are 


obtained for all (T;,A;) has already been discussed 


vo 








" oe 
7 a 
{smez3xe oft ovip doidw eovuis (Ajq onT .8L. £ ptt *~ 
be A, = 
Ss sgryetas entpszovs an ie esulsv . > 
Paver |" a 49 
Adsa) thet be B10VS ‘ens sone : a 
_sieb, 102 ae Ad ‘ean mi34648 - ma! 
« ’ 





Re RMN Age i a: 4 oe veer (i ae 
. ae ; ; 7 = . 


68 


maximum 


Fig. 3.19. The band of permissible models, in the 
p-A plane, for (T,A) data. The maximising 


p(A)<curve:sis: asp-shown:. 
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at the beginning of this chapter. By the same argument 
asetor 9(p,!,A) sets,) p= Mm; FOG A. eS Asvyd aye le 
maximise the integral, 7, and hence minimise the 
madvus sg). | At A(g) “there are two possibilities 


(1) A, < Ala) < 454) 


(2) hussein Sule 


t) 


Case "@S) Era’*ac? 4 Mist rated *in-tig 3 19 where; 'to*satisty 


travel-time conditions, 


oo, + u. Ae *whs:: 
ae) ee) eX 
Aver = J_1___ saul 
J J 


and Ds is chosen to maximise the integral - either at 
a e = 0 (active maximum) or at one of the 
Constraints on es (constrained or passive maximum). 

For the second case two local maxima exist which 
must be numerically tested to determine the global 
maximum. One is as described in (1) while the other 


is with p = : ae ele se interval and the same as (1) 


m 
in the Geunyee interval.= (0) 1s not essentially 
different from the general interval, except when uy =. 
As taking p = m,; will maximise il so per- 
turbing p as much as possible from the mean, while satis- 


fying travel time constraints will minimise ieee ene 
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However it is not possible to set 


Us A. < A. < 05 


(where ve is determined from the constraints) as at Asai 
there will exist a discontinuous increase in p from Wi 


to. u From the restriction dp/dA < 0 this increase 


ul cals 


is not possible and so the minimum must be constrained 


so that 
Py ay SS ihy Oo. 
p Ss 
Pitl Sipe ena itl 
where 
Rte Ty y PiAy z Pia Siel 
a ‘. A A bd 
see Sele al 
The values of the variables Pip i= Por. of, 7A UU See De 


determined numerically so that either a1,/8P; = 0 
(active minimum), p = u, (passive minimum with a1, /ap, <0) 
Or Pp = Was, (passive minimum with aI, /8p; > 0). This is 
illustrated an'fig. 3.20. 

At A(q)}, the contribution from the oe interval, 
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In determining the extremal average radii, the 
same arguments may be used as for (p,T,A) data sets 
to obtain the jp(A), which will maximise $¢(p,a,b)A(p) 
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where, as before m, 0° 5. 
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satisfy travel time constraints and maximise O( Diya, by (py). 


Thus 
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For Wi Fl shes mi, M, = m+ and Oo. = Asad while for 
m; < ish 4s Us, there are two possibilities:- m; =b 
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Fig. 3.21. The p(A) curve which will maximise the 
integral is shown as a solid line through 
the band structure, with b-a the interval 


over which an average radius is taken. 
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or ms > b, which must be tested numerically to obtain 
the absolute maximum of the iitedral. — A! so since in 
maximising the integral the Object is to inerease the 
area around p = b, taking ee = b will increase the 
value of @(p,a,b)A(p), as is shown in LiGae se cl. 

The values of p which will give the minimum 
d(p,a,b)A(p) curve are as shown in FAG. 32 22F where 


to satisfy time constraints, 


A * 

oT, BAL pF mPa 
“A * 
Eee so 


* A 9 0 4 , 
and m+ and Pj4 , are chosen to minimise the integral. 


Ingthe limit 


CL b = u; 4 Der 0. = Atal 
CZ) b = Usaye m= We, 05 = 05 : 


Obviously to program for this type f data is@more dif— 
ficult than, for (p,T,A) data so that though avprogram 


has been written to determine cp aeeenl (el and s (q) as a 


min 
LUNCtiOnN Of velocity fer (1,4) cate subroutines tor cind 
non-extremal models and calculate Rnax (212) and Rain (21?) 


have not yet been included. 
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the minimising p(A) curve, form (T,A) 


data for an averaging interval b-a, 
is represented by the solid line through 


the band structure. 
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ihe computatrons sacne thus far for’ (p,7T,4) have 
“useaq ideal) data in order to i2lustrate certain points 
about the general method of inversion. However once 
the analytical results of the following two chapters, on 
high and low velocity gradients with discrete data, have 
been programmed it will be possible to use real travel- 
time data to obtain an envelope of velocity-depth pro- 


files over the whole earth. 
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CHAPTER 4 


LOW VELOCITY CHANNELS 


AS was shown in Chapter 2, for both discrete and 
continuous data, it may be possible to distinguish between 
different types of discontinuity in the travel-time 
curve using amplitudes. For the type I or I' discontin- 
uity, the problem may be reduced by an ‘earth=stripping' 
process to a special case of discrete data and the con- 
ditions of maximisation and minimisation of velocity- 
depth models discussed in the previous chapter will be 
valid. Fig. 4.1 shows an example of such an earth- 
stripping procedure for the outer core where both P and 
S wave data were used (i.e. PKP-PcP, SKS-ScS) and the core 
radius was taken to be 3481 km. The diagram is merely 
intended to illustrate the type of results possible with 
ehnesprocedure.y Obviously if only PKP data is usedgthe 
extremal curves will not be as close. 


(k) 


However, for a type II discontinuity (i.e. nl(r +0) < 


Ae eo: IG Oh is undefined and the above method 
of reducing the problem cannot be used. Also it may not 
be possible, in practice, to distinguish between types of 
discontinuity from amplitude measurements, in which case 
the results of Gerver and Markeshevichs must be applied to 


discrete data. 
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McMechan (1971) has inverted discrete data con- 
taining low velocity zones, to obtain extremal velocity- 
depth profiles. However his method was to randomly 
Change the velocity function and depth of a low velocity 
channel to obtain those which gave a maximum or minimum 
velocity as a function of depth. It should be possible 
to obtain an analytic solution to this problem (i.e. 


(29) 


finding extremal integrals for q < p imvequatvone (2e 7000. 
This has been looked at by Chapman (personal communication) 
and will be discussed here, though no actual computations 


have yet been done. 


(1) 





For *q =% , the standard Herglotz-Wiechert inte- 
glamenay=vemused as *1in iChapter*s2S#or di oh, 
: ry {*) 
k = 
r)(q) = ya | 2 tant dr (4.1) 
k=1 
7 (kK) 


Where r‘*) is the top and r'K) ss the bottom of the k™ 


low velocity zone) must also be determined. This is 


obviously more difficult to solve directly than when 


= (k) (k) 


q 7 cur as Dotn FG are unknown. 


and 37 (3%) “tor © < © 


However, it is possible to place bounds on n(r) within 


(k) 


the low velocity channel, subject to the 6T and 6A 'K) (p) 


constraints. 
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To determine the range of velocities, it is 
necessary to determine the possible radii, s(q), where 


(by rewriting equation (2.7)) 


(o) k 
i 2a (q) i (k) 
s(q) = yay neal Se oc ae (a) | 


A low velocity channel may be constrained such that 


peer) <* 7) (0) as ali) ie for r a he Ve ie 


n - max 


where, by definition of a lcw velocity channel, need > ps 
(k) 


and ee sce aa for zero velocity. It can be shown (Appendix 
Ss) constraint, n(r) = ae) ee 


and ynl(r) = eran. for ie <r eas respectively maxi- 
(k) 


Ss)e that: using only the 67 


mise and minimise r = s(q) for most q < p 
It should be noted that in general,a model will 
only be physically valid if ds/dq 2 0, where the limits 
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as eucyogetend ito *sand Oucorrespond to q = p and a 


reflection, respectively. This implies that 
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which is negative, unless A(q) is on a reversed branch. 


However 








which is positive but not significant with respect to 


(kK) Thus to determine if a 


Ade CxXCePE LOR die. 9p 
towevefocrty channel is*physicalrly valid it is “only 


necessary to “check that 








dq 3 dq 


In’ the “case where a high velocity zone is immediately 
below the low velocity channel, az '°) sag will be small 
Gor’ trays SOctoming “fn tliis "regio and “So at ©) sag will 


be severely restricted from above. 


In determining the actual values of ieee and 
ie some other geophysical data must be used. Knowing 
ete Jay pe and piialh: the extremal values of the radius 
min ma 


at the bottom of the channel, are defined and for q < wat) 


Bie Ve tocity Woullas "are defined Dy i — eee and 


(kx) 


gut lied wad nc However, below a certain depth these low 


if) == 
velocity channel models will not define the extremal 
profiles and to do so the channel model must be slightly 


perturbed (see Appendix 3). It can be shown that 
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fractional variation in permitted velocities below a 
low velocity zone decreases away from the channel (i.e. 
the effect of the channel is most important in the 


vicinity of the channel) and 


(k) (k) 
Enax <a! yal) 
ti) (Gg) t'*) (9 OE? Saag AaEe 
min q Ase ) 


So far the sa") (p) constraint has not been fully 
utilised as, for example, if nt®) (x) + o, the a oh? dis- 
continuity may be satisfied by an infinitesimally thin 
low velocity channel but this would give sa ‘*) (p) + 0 


which would not agree with the actual data. Also gx 


can be satisfied with a thick channel and n(r) > gel. 
but this would create a large jump in sa tp), which 
again is incompatible with the data. Therefore it should 
be possible to further restrict low velocity models using 
the 50%) (p) constraint. 

As stated in Chapter 2, low velocity models may 
have numerical instabilities and to avoid this the turn- 
ing radius may be averaged over a range of n(r), as shown 
by Bessonova et al, 1970. Maximising and minimising an 
averaged radius for discrete data, without low velocity 


zones, was discussed in the previous chapter. The inclu- 


sion of low velocity channels will be briefly discussed 
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Neren Slt wile be recalled etnat the cerinition Of 4 radius 


averaged over some interval a < p < b is 


n b 
R(a,b) 1 2 k 
co arene noe (| A (P) 4 (p,a,b) dp + | ee 1) (q)dq 
o k=1 
a a 
(Ar 2) 
where Toe a) is as defined in equation (4.1) and 
-la -l b 
Cc oe Drs —_ b < 
OSs p COs Dp ee) ne 
o(p,a,b) = 
-la 
= a< <b : 
cos p p 
As 
ene He 
ope ae 
k 
PMs as A(p)dp + y¢ or (4.3) 
k=1 
No 
and hence, integrating by parts 
Oo O No 
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a a a 
n 
ome) 
ap 
+ | | A (q) dq 3p dp 
a p 
i k 
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using equation (4.3) and the condition ¢(a,a,b) = 0. 


Thus equation (4.3) may be rewritten as 


ia} 
Riavbip ee 1 3 eer (k k 
a ae ss | aD 3p ap+ p= JF oro ip") ,a,b) 
k= 
5 a 
b 
k 
is 
= sts | yo 8) (aq) aq (4.4) 
k=1 
a 
where 
a b 
sop agp leh <5 els ii 
2 i a 
B(paaas a tp-ebo) = i 
cy = 
op 


Finding the extremal values of R(a,b) then becomes_a 


question of bounds on Tt. Because the sign of ¢* (p,ayb) 


varies across p = b, T = Trin \P? will maximise radius for 

qo- pea eb Wwitle LOmepDe< pi< Nor T= tn eed e will maximise 

the radius. “~The values of 7 and t_.. will depend on 
max min 


the data set. For the simple case of T; T(A;) and 
Sane p(A.), the upper and lower limits of T, in the wip 
plane, become a box-structure as shown in fig. 4.2. 


Obviously in this case 
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ray parameter p, sec/deg 


The upper and lower limits on 7 for" (T,p) data, 


shown in the t-p plane. 
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and given other data sets, T may be suitably limited. 
The second term on the right hand side of equa- 


Bee: 


tion (4.4) may be estimated, with limits, if (T., p. 4 


alt 


are measured. At worst, it must be assumed that there 


are discontinuities in every interval and then 


(ki) 3 m4 
6t et: MES 
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with p such that 9(p ,a,ob) is maximised in that 
interval. 

The third term on the right hand side of equation 
(4.4) may also be maximised and minimised using T(p) and 


gr (kK) data, 
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CHAPTER > 


TRIPLICATIONS 


A triplication on the travel-time curve is caused 
Dyan increase ian veloGity gradient. “For a triplication, 
themconditiion dp/dA = 0risno longer valid, as. for the 
reflected part of the T-A curve p becomes a monotonically 
increasing function of A. Thus (as can be seen in 
fig. 5.1) that segment of the Herglotz-Wiechert integral 
is negative. For continuous data this presents no pro- 
blem. However with discrete data a number of possibilities 
must be considered. 

If triplications are assumed to exist, then: 

(1) the triplications occur between data points and are 
unobservedquaie) turiplications occur at data points but 
arewunobserved or (3) triplications are observed at) the 
data sports. 

The first case will not cause any change in the 
band structure. For the second case, the data point 
observed may be on the reflected branch such that Py is 
less than Psa] where the abs measurement is again on 
the forward branch. This situation is indicated in 
Eigen 2s. lethough of course the p (A) curve shown in this 
figure is for continuous data. If this is the case at A; 


then using the standard definitions of Mi, Uz and w, as 
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a triplication with an example of possible 
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points, superimposed. 
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in Chapter 3 (1.6. T;Atdata is being considered) then 
to allow for this possibility t:she band between A; and 
A;4, Must be widened (as illustrated in fig. 5.2). 

For the case of an observed triplication, which 
has been studied by McMechan (1971) and McMechan and 
Wiggins (1972) the time constraint can no longer be 
satisfied on certain branches, when maximising and 
minimising the integration curve. To maximise in the 
interval containing the triplication, therefore, it is 
possible to take p = Us, LOG, es u, (when again consi- 
dering the general case of T,A observations with no p 
data). However, compensation must be made for this 
branch's lack of agreement with travel time constraints 
at some stage. 

The maximum curve will be as shown in fig. 5.4 
where the correspondence between the branches of the 
T-A curve and the p(A) curve at a triplication are as 
shown in figs. 5.3 and 5.4. The section labelled A in 
both diagrams refersto the forward branch, as does that 
labelled C. The section labelled B refers to the reverse 
branch of the travel time curve. The reflected branch 
is fixed either by time constraints or such that it 
passes through (A;, W,)- The point p', shown in the 
diagram, defines the value of q for which this curve is 


an absolute maximum. More generally, the maximum for any 
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The band structure, in the p-A plane, for an 
unobserved triplication at Aj. Thescross 
hatched area shows the usual band width if 

no triplications occur while the shaded region 
shows the increase in band width necessary to 
accommodate the type of situation illustrated 
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Portion of the T-A which includes a 
triplication. The labels on the 
difference branches of the curve 
correspond to sections of the p(A) 


curve shown in fig. 5.4. 
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Thesp (A) scurve in the Ve interval which 
will maximise the Herglotz-Wiechert integral 
when a triplication is observed at A.- 
Labels A and C refer to the forward branches 
of the travel-time curve, as shown in fig. 


5.3 while B refers to the reverse branch. 
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particular gq may be obtained from fig. 5.5, where 2) is 
such that travel-time constraints are satisfied and the 
integral TG) is a maximum. Below a certain value 
Ofed, though gq is still greater tian Wir the curve of 
integration is as shown ingfig. 5.6. gin the limit as 
ae Wr the reverse branch disappears, Da > ™m, and the 
case becomes that discussed in Chapter 3. 

Theiminimisation [omqe- mM; igmas shown in fig: 
be. / jawoere De =q. For-q qm, it is not obvious whether 
fig. 5.8(a) or 5.8(b) would provide the absolute minimum, 
though this may be tested numerically. AS q7W,, the 
case is again that discussed in Chapter 3, for an absolute 
minimum. 

When triplications are observed at data points 
it is possible to include the maximising and minimising 
procedure necessary at a triplication in the more general 
computations considered in Chapter 3. However, this has 
not been done in the computations of this thesis as it 


has already been discussed by McMechan (1971). 
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The more general maximisation curve 


in the jth 


interval, where DP; is the 
value of p which will maximise t ©) (aq) 


and also satisfy time constraints. 
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Fig: 


5.0 ()i. 


Both this and fig. 5.8(a) show 


possible minimisation curves for 
q< m,. The curve giving the 
absolute minimum may be determined 


numerically. 
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CHAPTER 6 


CONCLUSIONS 


The results of the experiments with sampling 
intervals show that the more data available, the 
narrower the limits on the acceptable velocity-depth 
models become, which is as expected. piual ds by 
fixing part of the model, it appears that though 
making more assumptions about the higher regions of 
the earth does further restrict possible models in 
the near vicinity of the known Pod iGat Tie change 
in envelope width occurs at: lower depths. 

The results of computing individual profiles 
and their behaviour with respect to the extremal 
curves shows that although these profiles do touch 
the extremal curves at various depths, in general they 
lie within the envelope (though their exact behaviour 
depends on the actual data). In fact, the extremal 
curves,which show the largest ambiguity in velocity 
possible at any particular radius, are not themselves 
physically valid profiles. Thus averaging these extremal 
curves is not meaningful and to see the effect which 
averaging would have on the permissible amount of velo- 
city variation, an averaged radius is defined and then 
suitably maximised and minimised for various averaging 


lengths. 
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Computations of these averaged envelopes for 
a number of averaging lengths show that averaging 
does reduce the envelope width considerably and for 
the range of depths considered the envelope width is 
reduced by 40% for an averaging length of 0.4 sec/deg. 
The exact way in which extremal curves vary with 
averaging length is dependent on the actual data 
(for similar reasons as for the exact behaviour of 
particular profiles with respect to the extremal curves, 
mentioned above). This type of computation is essential 
for any comparison with the more indirect inversion 
methods, such as Gilbert and Backus or Monte Carlo. 

Maximisation and minimisation when high or low 
velocity zones are included is discussed in this thesis, 
though the analytic results obtained here have not been 
programmed. When this has been done, it will be possible 
to invert real data, which will facilitate a qualitative 
comparison of the extremal averaged curves obtained by 
this inversion method with the results of other inver- 
sions. For an interface such as the core-mantle boundary, 
it. is possible to remove the effect of the low velocity 
zone by an ‘'earth-stripping' technique which effectively 
removes the earth's mantle and then inverts the reduced 
data. This technique is illustrated in this work, using 


both p and s wave data. 
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So far, to compute one velocity-depth envelope 
takes about 30 secs C.P.U. time, though when calcula- 
tions for high and low velocity zones have been 
inemided, this time will sobviously sincrease,., (A 
program for (T,A) data exists but does not contain 
subroutines to determine extremal averaged curves. 

The necessary programming for (p,A) data would be 
fairly simple and it is expected that both (T,A) and 
(p,A) envelopes will be wider than a corresponding 
(Dio eoenve lope. 

With the present program, it should be possible 
tomstudyethe efrect on the profiles of errors inthe 
measurements, though a thorough analysis would require 
a systematic variation in all the given data parameters. 

One of the more important aspects of this thesis 
is that the very nature of discrete data implies ambi- 
guities in any inversion procedure. Tests show that 
these ambiguities may be recluced by increasing the 
amount of data and it may be further deduced that using 
time, distance and slowness data simultaneously will 
give the most restrictive envelope. However, there will 
always be some ambiguity. 

The other important con¢lusion reached is that 
it is essential, when comparing with results from other 


inversion techniques, to understand exactly how these 
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comparisons should be made. The concept of an averaged 
radius, which must be maximised and minimised, is in- 
troduced so that it will now be possible to compare 
results from a Herglotz-Wiechert inversion procedure 


with other types of inversion. 
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APPENDIX 1 


At high frequencies (Richards, 1971) the wave function 


¢ is described by the Helmholtz wave equation 
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Using the above results, Chapman (1972) has shown the 


following: 


A iip tice discontinuity: 
(k) 


Orsi. ; 5? (r) may be approximated by a 
linear function, so that equation (Al.1) reduces to 


Stokes equation. Assuming the free boundary condition, 


6 (r ) = 0, the reflection coefficient becomes 
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R(r_,v) (2D + dt (p)) Mies to! 
PUNE eee Lat AB) ie 17 3 (Al.2) 
- P pantera s/ 7) 
where 
6k) 
ca t2/3 = -i | Er) dri. 
S (p) 


((> © at no interface present). 


The poles of R(ryrv) exist at 
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for a diffracted signal with A > A , where 
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It should be noted that at large ranges, the j = 0 term 
will predominate where ee 2.3381. More realistic 
boundary conditions cause poles intermediate between 


the Airy zeros but with the same frequency dependence. 


Type II discontinuity: 
nite ie?) = 0 


Therefore, for r > p(k) E* (r) must be approximated by 
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a quadratic. Equation (Al.1) then has the parabolic 
cylindrical functions as solutions. The reflection 


coefficient is (Phinney, 1969) 


_e ee (a,—t) ; 
R(ro,V) = a CP exp{iwt (p) } (Al. 4) 
where 
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which is independent of frequency. Again, at large ranges 


only the j = 0 term will be important. 
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R(xro7¥) will have poles where 
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be asymmetric. The exact numerical factor depends on 
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APPENDIX 2 


For a™typevliediiscontinuity (Chapman, 1972)m@the 
range of the ray tends to infinity as p tends to oot). 
That is beyond a certain A(p), the original becomes 
infinitely small in amplitude due to increased geome- 


trical spreading. Consider 
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To determine the amplitude, which is related to 
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and hence the amplitude is governed by 
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This result is identical to the asymptotic form of the 

wave solution (Al1.5 with j 
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APPENDIX 3 


As shown by Chapman (1972), n(r) = Te and 
k . 
n(r) = gees will respectively maximise and minimise 


r = s(q). From Chapter 2, it will be recalled that 
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To determine whether (A3.3) gives the absolute minimum 
integral, n(r) is perturbed to Mew (6) + 6n(r) and the 
corresponding perturbations in z (x) and t'*) (q) observed. 
By differentiating equation (A3.2) and substituting from 
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Obviously, to determine whether n(r) = ieee maximises 


the integral it is necessary to determine the sign of 
the bracketed terms in equation (A3.5), and this may 
be done on a qualitative basis by plotting the behaviour 
of these terms. -biei-t-—tse-first—assumed that ren has a 
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channel, then 
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These approximations to the terms in brackets in equation 
(A3.5) are as shown in fig. (A3.1) where it can be seen 
that for the most part, the term in brackets is negative, 
making sr) ¢q) negative for positive On. That is, that 
n(r) = n &*) (x) will maximise es ) qe and, therefore, 
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minimise velocity for any given q. The question of 


maximising velocity may be similarly treated, with n(r) = 
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are then perturbed by changing n(r) to qigeees! - 6n(r). 
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However it should be noted that at the points 
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a and b shown in fig. A3.1 the two curves inter- 


sect so that the bracket term in equation (A3.5) changes 
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The exact perturbation of the low velocity channel which 

is required is not obvious and should probably be deter- 

mined by trial and error. This will not be too ambiguous 

because both the q and r ranges are restricting factors. 
The above discussion is probably best illustrated 

in the radius-velocity plot of fig. A3.2, where it should 

be noted that in a region below the low velocity channel 

the channel model does not correspond to the extremal 

profiles. These must be obtained by perturbing the channel 


model within the shaded area of the diagram. 
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C3) 





Piggihoee. A radius-velocity plot for a low velocity 
channel with ee and eee the maximum and 
max min 
minimum radii for the bottom of the low 
velocity channel (as determined by gy 
The shaded curves are those obtained from 
the low velocity channel model while the 


cross-hatched area is produced by perturbing 


the model. 
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